A model differential equation for turbulence 
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A phenomenological turbulence model in which the energy spectrum obeys a nonlinear diffusion 
equation is presented. This equation respects the scaling properties of the original Navier-Stokes 
equations and it has the Kolmogorov —5/3 cascade and the thermodynamic equilibrium spectra as 
exact steady state solutions. The general steady state in this model contains a nonlinear mixture 
of the constant-flux and thermodynamic components. Such "warm cascade" solutions describe the 
bottleneck phenomenon of spectrum stagnation near the dissipative scale. Self-similar solutions 
describing a finite-time formation of steady cascades are analysed and found to exhibit nontrivial 
scaling behaviour. 

PACS numbers: 47.27.Eq, 47.27.Gs 

I. MODEL EQUATION 

In this letter, we present a model equation for the energy spectrum of isotropic homogeneous turbulence, 

Et = i(A:"/2^i/2(i?/fc2)fe), + / - ^k'E, (1) 

where t is time, k is the absolute value of the wavenumber, v is the kinematic viscosity coefRcient, / is an external 
forcing and the energy spectrum E{k, t) is normalised so that the kinetic energy density is / E dk. The third term on 
the RHS of this equation is obvious and describes the linear process of viscous dissipation. The first term contains a 
nontrivial model of the nonlinear processes in turbulence which rests on three basic properties 

1. The total energy density J E dk is conserved by the inviscid dynamics. The characteristic time of the spectral 
energy redistribution for local interaction of turbulent scales is of order of the vortex turnover time, l/Vk^E. 

2. The steady state in forced turbulence corresponds to a constant energy cascade though the inertial range of 
scales which is described by the Kolmogorov spectrum, 

E^CP^'^k-^l^, (2) 

where P is the energy flux (constant in t and k) and C is the Kolmogorov constant. Experimental measurements 
^ : give C = 1.6 ± 0.17. As we will see below, equation (QJ has an exact solution of form (0) with C = (24/11)^/^ « 1.68 
^ (C can be changed by tuning the constant factor in the first term). 

, 3. When the wave- number range is truncated to a finite value and both forcing and dissipation are absent, turbulence 
reaches a thermodynamic equilibrium state characterized by equipartition of energy over the wave- number space 
O ^' In terms of the one-dimensional energy spectrum this means E (x k'^ which is obviously a steady state solution of 
the equation for f = v = 0. Mathematical simplicity and respect to the above basic properties of Navier-Stokes 
turbulence make model Q useful for practical and numerical applications. We will now analyze solutions of in 
greater detail to find other properties that are predicted by this model. 
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II. STATIONARY SOLUTIONS 



Let us consider steady-state spectra in the inertial range. For f — v — 0, we have the following general time- 
independent solution, 

i? = CA:2(Ffc-"/2+Q)2/3, (3) 

where C — (24/11)^/"^ « 1.68 and P and Q are arbitrary constants. For Q = 0, this gives the pure Kolmogorov 
cascade solution , whereas for P = this is a pure thermodynamic equilibrium. For the general solution, both the 
constant flux of energy P = — ifc-'^^/^i?^/^(£'/fc^)fc ^ and a thermodynamic part (Q 7^ 0) are present; they appear as 
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an interesting nonlinear combination and not just as a linear superposition because equation |^ is nonlinear. Thus, 
one can refer to solution ^ with finite P and Q as a warm cascade to distinguish it from the pure Kolmogorov 
solution which could be viewed as a cold cascade. 

Let us suppose that turbulence is produced near some scale k — kg so that f{k) > only in a finite range in the 
vicinity of fco- Suppose that f{k) — to the left of this range (at large scales) and in a large inertial range to the 
right which ends at a very high k ^ kd where viscosity or some other dissipation mechanisms (/(fc) < 0) become 
important. Then, up-scale of the forcing there will be a pure thermodynamic solution with P = and Q because 
there is no dissipation or forcing assumed to be present near fc = to absorb or generate a finite energy flux. In 
the inertial range down-scale of the forcing there will be a constant flux cascade solution. This solution typically 
takes the form of a pure Kolmogorov (cold) cascade and extends down to the dissipation range where the energy 
flux is absorbed. Typically, the solution only penetrates a finite distance into the dissipation range and adapts itself 
until it provides sufhcient dissipation to absorb the supplied flux. In the presence of dissipation, the model does not 
develop structure at arbitrarily high k as it would in the inviscid case. The qualitative features of the steady state 
are independent of the detailed form chosen for the dissipation. Figure d shows the steady state solutions obtained 
numerically for several different choices of the dissipation. 

However, if the dissipation is not sufRciently strong, the solution can penetrate far enough into the dissipation 
range to reach the maximal wave-number which necessarily exists in any numerical solution. If one imposes a zero 
flux condition at the right end of the computational interval, the energy flux is reflected from the maximal wave-number 
leading to greater values of E in the dissipative range. Such a cascade stagnation acts to enhance the dissipation rate 
and thereby to adjust it to the energy flux to be absorbed. This phenomenon is common in numerical simulations of 
turbulence and is usually called the bottleneck phenomenon 3]. Figure |21 shows a numerically obtained steady state 
for the dissipation function. 



with 1^0 ~ 1.0 X 10^^ and kjj — 500. The bottleneck phenomenon is clearly seen as an energy "pile up" over the 
cold cascade solution near the dissipative scale. In our model, the bottleneck phenomenon is described by the warm 
cascade solutions; in particular the theoretical curve in figure|21is computed by taking P « 14.5, Q « 1.5 x 10^^ in 
equation The relative importance of the "thermal" effects with respect to the cascade grows as one moves from 
larger to smaller scales; in extreme cases the spectrum can be nearly pure Kolmogorov near the forcing range and 
almost purely thermodynamic near the dissipative scale. 

Because of the nonlinearity, solutions for given forms of forcing and dissipation are usually hard to find analytically 
and one has to use numerics. However, a lot of insight about the qualitative behavior of the system can be gained 
from considering stationary solutions Q in an inertial range ki < k < k2 and fixing the spectrum at its boundaries, 
E{ki) = Ei,E{k2) = E2. This kind of the boundary conditions roughly models the forcing and the dissipation effects 
outside of the inertial range. This gives 



Thus, the sign of P is opposite to the sign of {E2/E1 — fc^/fc^) and can be either positive or negative depending of the 
spectrum steepness with thermodynamic k^ solution been a borderline case for which P = 0. Constant Q can also be 
either positive or negative with Kolmogorov —5/3 been the borderline slope. It is convenient to think of the Q < 
solutions as negative temperature states. 



So far we concentrated on the stationary solutions but how do these solutions form? We consider from now on the 
inviscid case. The Kolmogorov —5/3 spectrum is of a finite capacity type in that it contains only a finite amount 
of energy at the high k end. Let us take an initial condition which is compactly supported and force the system by 
imposing a constant fiux boundary condition across the left end of the computational interval. Owing to the finite 
capacity of the Kolmogorov solution, an infinitely remote dissipative scale must be reached in a finite time. The 
solution has a nonlinear front at A; = k^{t) and this front accelerates explosively, reaching fc = 00 at a finite singular 
time which we shall denote by t^. We can equally consider the decaying case where an initial distribution of energy, 
compactly supported at large scales, is allowed to spread under the action of the nonlinearity without any external 
forcing. Provided there is sufficient energy in the initial condition, the right front still reaches k = 00 within a finite 



z/(fc) — voik — fcu)^ k > kjj 
fc < fc_D 



p = m/ckif/' {E,/ckif'^]/{k-''^' - fcr"/^) 

Q ^ [kl'^(E2/Cf" - kl'-\E,/Cf'^]/{kl"^ - k\"^) 



(4) 
(5) 



III. NON-STATIONARY SOLUTIONS 
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steady states obtained for different dissipation functions 




FIG. 1: Numerically computed steady states for several choices of dissipation function: iyi{k) = (500 — k)^ , V2{k) — 4.0 x 
10"''(500 - kf, va{k) = 1.0 x 10"^fc^ The Kolmo gorov spectrum is also shown for comparison but shifted slightly for clarity. 



Evoiution of bottieneck for v{k) = 1 E-5 (k-500)''2 




FIG. 2: Numerical evolution of a bottleneck for dissipation function !/(fc) = 1.0 x 10 ^{k — 500)^. The resulting steady state is 
well described by the solution ^ . 



time but in addition a second front propagates to the left spreading the spectrum to large scales. This second front 
does not exhibit singular behaviour in finite time. Hence the decaying case is broadly similar to the forced case as far 
as large k behaviour is concerned since the initial concentration of energy at large scales acts as an effective forcing 
for the right front. 

It is well known (see for example ,4J) that the solutions of nonlinear diffusion equations with compactly supported 
initial data often have the property of remaining compactly supported during the time evolution. This turns out 
to be the case here. The left and right nonlinear fronts actually correspond to the left and right extrema of the 
support of the solution. Since these points must be determined as part of the solution, we are, in principle, required 
to solve a moving boundary problem with two free boundaries. In order that the problem be well-posed, we require 
an additional movinq boundary condmon on each interface's?. For the problem under consideration, the appropriate 
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FIG. 3: Integration contour leading to the moving boundary condition, ||SJ. 



condition for the right front is 



dk 



(R) 



dt 



hm k^E~i— (k-^E) 
(«) dk ^ ' 



(6) 



A similar condition holds for the left front. 

This formula can be obtained as follows. Consider the time interval [0,t] with t < t* . For < r < t the right 

interface is given by the curve ki^\T) as shown in figure El Taking, e << 1 and assuming that the solution goes 
continuously to zero at the right interface we define the curve k^ir) by the condition E{k^{T),T) = e for < t < i. 
We now integrate around the contour, T, in (k, t) plane as shown in the figure, 



We obtain 



Proceed as follows 



dE dP 
df'^'dk 



dk dr = 0. 



(j) (Edk + Pdr) = 0. 



(7) 



(8) 



E (dk + Pdr) = 0, 



(9) 



where 



P^-k'^E--^-^(k-'E) 
ok ^ ' 



(10) 



fc<«'(t) 



Edk 



E dk + Pdr 



/Cc(O) 



Edk 



Now E' = on ki'^^ (r) and E = e on ke (t) . Using these facts we obtain 



fe,(r) 



E {dk + Pdr^ ^ 0. 



h{t)-kM 



PdT= I —dk- 



fe.(t) 



-dk. 

e 



(11) 



(12) 
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We now take e ^ so that fce(r) ki^^r). The integrands on the RHS are bounded by 1 and therefore give no 
contribution. We are left with 

ki"^^ {t) = /fc(^) (0) - /* Piki"^^ (r) , r) dT, (13) 



^0 

which in conjunction with pU|l yields (jSJ upon differentiation. 

Let us now look for a self similar solution of equation taking the following form 

E = {U-tYF{r^)- r] = k/k,, k^ ^ c{U - tf , (14) 

where a, h and c are constants. Clearly, h must be negative since we require that fc, — > oo as i — > t,. Substituting p4|) 
into with f = v = Q wc find that the t dependence drops out of the equation if a = —2 — 36. We then have the 
following equation for F , 

(-13/2 

(36 + 2)F + hr^F' = r_(r;"/2^i/2(^/^2yy^ (^5^ 
8 

where prime means differentiation with respect to 77. Equation (|15|l defines a one-parameter family of self-similar 
solutions. Note that substitution of the form (|14|l into the moving boundary condition, lO yields the same similarity 
relations. Thus the moving boundary condition is consistent with the similarity ansatz but does not provide any 
additional constraints. In particular it does not tell us which member of the 1 parameter family of solutions is 
selected by the PDE. 

The solution near the front tip can be found by expanding F in series with respect to small (1 — rj); in the leading 
order we have 

1662 

F=-^{l~ilf (16) 

which gives for the spectrum 

1662 / fc ^ ^ 



kJ'iU-t)^ \ kj 

Notice that the quadratic decay of the solution as one approaches the tip is exactly the decay required in order to 
ensure that the speed of the tip as given by (j^ is finite and non-zero. 

We are interested in solutions which behave like a power law far behind the front. That is E k^^ as fc ^ 0. 
The relations H14|l then imply that x — —a/b. The pure Kolmogorov spectrum, x = 5/3, therefore requires 6 — —3/2, 
corresponding to what one might consider to be normal scaling in the wake of the front. 

We performed numerical simulations of the forced solutions of equation JQ) with compact initial data and constant 
input flux on the left. We use the numerical solution to check the development of a self-similar front with tip of form 
(|16|l and to determine which value of 6 is selected. The results are shown in figure figure^ The scaling parameter, 6 
and the singular time, t^, are most conveniently extracted from the relation 

which allows one to calculate 6 and from a linear fit of the data near t = as shown in figure |31 We find that 
= 0.0799 and 6 — —1.748 which corresponds a significantly steeper than Kolmogorov slope, x — 1.856. The singular 
time, depends on the choice of initial conditions but the anomalous scaling exponent does not. In particular we 
verified that the same value of 6 is obtained for the decaying case. Such anomalous scaling behavior whereby the 
exponent of the solution in the wake of the nonlinear front is larger than the Kolmogorov value has been observed 
before. Examples include MHD wave turbulence 6] and optical turbulence 7]. We expect that this anomaly is a 
general property of finite capacity systems rather than a property of our model and therefore may also be present in 
the Navier-Stokes equations. 



IV. ORIGIN OF THE TRANSIENT SPECTRUM 



For the model the origin of the anomaly can be traced to the question of existence of a solution of the similarity 
equation H15|l which has the correct behaviour both for the wake and the tip. Written in terms of x rather than 6 and 
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FIG. 4: Forced time-dependent solutions beginning from compact initial data showing development of self-similar front with 
power law wake. 
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FIG. 5: Calculation of the asymptotic scaling properties of the self-similar solution. 



re-scaled to get rid of constants, the similarity equation is 

2 f dF \ d / 11 ^ d , , A 

We require that this equation have a solution which behaves as 77"^ as 77 ^ and also satisfies the boundary condition, 
©, at the front tip, 77 ^ 1. Such a solution is not typical and actually exists only for one value of x. In particular, 
such a solution does not exist for x = 5/3, the Kolmogorov value. The structure of the problem can be studied by 
introducing a new independent variable, s = log 77 and a pair of dependent variables, /(s), g{s) defined by 



(20) 
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Trajectories in {Kh) Plane 
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x=1.86 
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Fixed points 



FIG. 6: Trajectories in the {k, h) plane. 



The purpose of this transformation is to "autonomisc" the equation 10] . That is, we remove the exphcit dependence 
on 77 from the equation. When this is done, equation l|19() is equivalent to the following autonomous first order system: 



ds 



(21) 
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5/' + 6/.g- V + -{if + xg) 

X — 6 



and methods of phase plane analysis can be applied. The associated dynamical system has three equilibria PI = (0, 0), 
P2 = (0, 10/3(a; — 3)) and P3 = (1, —1). Note that PI and P2 are singular points of the original equations, H21|) . 

Let us now attempt to identify the trajectory in the (/, g) plane describing the physical solution. We note that 
the required trajectory must remain in the quadrant / > 0, g < 0. The wake is at 77 —s- 0, or s —s- 00. In the wake, 
F ^ rj^^ and F' ^ —xri~^~^. It follows that / ~ 77^"^ and g ~ —x'q'^~^. The value of x is less than 3 so both / and 
g go to zero as 77 ^ 0. The wake is therefore at (0, and should be reached as s — > —oo. 

Now consider the tip. We know that the front tip is at 77 = 1, or s = O.Both F and F' go to zero as 77 ^ 1 but 
it follows from © that the product F~^^^F' remains finite. The front tip therefore lies at some finite point on the 
negative g axis. In fact the tip lies at P2. To show this, we show that the moving boundary condition is satisfied 
there. Translating (jSJ into the self-similar variables we obtain 



hT' 



b-l 



= h+ia 1 I dF 

-r2^2 lim 772^— — . 
■n^i dr] 



(22) 



Balancing the t dependence gives us back the self-similarity condition, a 
changes of variables and rescalings to obtain 

3 

h = lim —h. 
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-2 — 36. We now carefully do the various 



(23) 



Noting the similarity relations (|14|) . we see that this condition is satisfied identically at the point P2- 

The critical value of x occurs when the unstable manifold of Pi intersects the stable manifold of P2 thus forming 
a connection between the two singular points of the original equations. This is then the only trajectory which can 
satisfy the required conditions to describe both the tip and the wake. A numerical approximation to the unstable 
manifold of the origin as x is varied through the critical value is shown in figure |^ The corresponding front profiles, 
converted back into the self-similar variables are shown in figured 

In practice, any model should include dissipation so that this self-similar solution above will be valid only until 
the front tip meets the dissipation scale. After this, the transient slope gets replaced in the inertial range by the 
stationary cascade solution, with or without bottleneck depending on the dissipation, as discussed above. 
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FIG. 7: Corresponding front profiles. 

In summary, we presented a simple model in which basic properties of the Navier-Stokes turbulence are built-in: 
dimensionality and scaling, Kolmogorov and thermodynamic spectra. This model allows to obtain new predictions 
about the "warm" cascade states which are a nonlinear mixture of the cascade and thermodynamic solutions and 
which describe the bottleneck phenomenon. The model also allowed us to study the self-similar dynamics of the 
finite-time formation of the steady cascade states. 

We thank Ildar Gabitov and Alan Newell for helpful discussions. 
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